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NONTECHNICAL  SUMMARY 


In  this  paper,  we  introduce  a  new  technique  for  analyzing 
simulations  of  stochastic  systems  in  the  steady  Btate.  From  the 
viewpoint  of  classical  statistics,  we  address  the  questions  of 
simulation  run  duration  and  of  starting  and  stopping  simulations. 

We  are  able  to  do  so  by  avoiding  two  difficulties  which  have 
previously  made  classical  statistics  Inappropriate  for  simulation 
analyses.  These  are  the  statistical  dependence  between  successive 
observations  and  the  inability  of  the  simulator  to  begin  the 
system  in  the  steady  state. 

For  many  stochastic  systems  being  simulated  it  is  possible 
to  find  a  random  grouping  of  observations  which  produces  Independent 
identically  distributed  blocks  from  the  start  of  the  simulation. 

This  grouping  then  enables  the  simulator  to  avoid  the  two  problems 
mentioned  above.  He  has  at  his  disposal  the  methods  of  classical 
statistics  such  as  confidence  intervals,  hypothesis  testing, 
regression,  and  sequential  estimation  which  are  appropriate  for 
independent  observations.  Furthermore,  information  that  is  useful 
in  estimating  the  steady-state  behavior  of  the  system  can  be 
collected  from  scratch  thus  eliminating  the  problem  of  the  initial 
transient . 

The  approach  mentioned  above  is  appropriate  for  the 
simulation  of  systems  which  returns  infinitely  often  to  a  single 
state.  In  the  current  paper,  we  restrict  our  discussion  to  the 
general  multi-server  queue,  with  arbitrarily  distributed  inter-arrival 
and  service  times.  We  leave  the  more  general  systems  to  future 
publications . 
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Section  2  reviews  some  relevant  results  from  the  theory  of 
statistical  inference.  Sections  3  and  4  discuss  the  probabilistic 
structure  of  stable  queueing  systems.  A  queueing  system  is  stable 
if  the  customer  arrival  rate  is  strictly  less  than  the  maximum 
service  rate  (with  all  servers  working).  In  this  case,  it  •'ay  be 
shown  under  mild  conditions  that  the  idle  state  (the  state  in  which 
all  servers  are  idle)  occurs  Infinitely  often.  Furthermore,  letting 
a  busy  cycle  refer  to  the  time  interval  between  two  successive 
idle  states,  it  may  be  shown  that  observations  made  in  different 
busy  cycles  are  statistically  independent  and  identically  distributed. 
These  observations  might  b?,  for  example,  the  number  of  customers 
served  in  the  busy  cycle,  the  length  of  the  buoy  cycle,  the  sum  of 
the  waiting  times  of  customers  served  in  the  busy  cycle,  or  the  sum 
of  some  function  of  the  waiting  timer  of  th«  customers  served  in 
the  busy  cycle. 

It  is  shown  that  each  busy  cycle  provides  information 
relating  to  the  system* 3  steady-otcte  behavior.  In  particular,  the 
expected  value  of  any  well-behaved  function  of  the  steady-state 
waiting  time  is  equal  to  the  expected  value  of  the  sum  of  that  same 
function  of  the  Individual  customer  waiting  times  in  a  single  busy 
cycle,  divided  by  the  expected  number  of  customers  served  in  a  busy 
cycle.  This  fact  may  be  uoed  with  the  above-mentioned  independence 
to  enable  the  simulator  to  perform  a  thorough  statistical  analysis 
of  the  steady  state.  One  merely  directs  his  analysis  toward  the 
estimation  of  properties  of  the  individual  busy  cycle  and  then 
Infers  corresponding  properties  for  the  steady  state.  The  former 
task  is  simplified  because  of  the  independence  and  identical 
probabilistic  structure  of  different  busy  cycles,  which  permits 
classical  statistical  analysis. 
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To  illustrate  the  application  of  these  ideas,  consider  the 

problem  of  obtaining  a  confidence  interval  for  the  expected  steady- 

state  waiting  time*  K{W}.  Let  denote  the  sum  of  the  customer 

2 

waiting  times  in  the  kth  busy  cycle  and  let  i  and  denote 

respectively  the  sample  mean  and  sample  variance  for  iri  N 

observations  (busy  cycles).  Let  denote  the  number  of  customers 

2 

served  in  the  kth  busy  cycle,  and  let  a  and  sq  denote 
respectively  the  sample  mean  and  sample  variance  for  in  N 

observations.  To  obtain  a  confidence  Interval  for  E{W}  with  at 
least  100(1-y)%  confidence,  we  obser  *e  the  system  for  N  busy 
cycles.  The  interval  computed  is  then 


max{0,  /2  SyAfi} 

"  +  Zl-y„/2 


E{W} 


5  +  *1-^/2  *Y//B 


'"“'O-  “-'l-v2/28a 


/*^} 


where  zx  is  the  lOOx  percentile  for  the  norv-V-.l  distribution  and 

(]L  and  y2  satisfy  Yj_  +  Y2  *  Y  - 

Section  5,  a  two- >  cep  procedure  is  developed  for  obtaining 
an  approximate  confidence  interval  for  the  -?teady-8>:ate  mean  of  a 
»?,enerai  function  of  the  customer  waiting  time.  The  first  stc\ 
consists  of  a  short  simulation  run  which  serves  as  a  planning  guide 
for  a  much  longer  seccjkv  tun.  4c  the  end  of  the  first  run,  the 
simulator  is  able  to  estimate  the  ntf-iimum  confidence  interval  length 
which  ran  be  obtained  from  the  second  run  c,iven  specified  level  of 
confidence  and  run  time.  The  procedure  is  thus  viewed  as  a  tool 
by  which  the  simulator  may.  balance  computation  cost  against  level 
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of  precision  and  hence  implement  a  rational  design  of  experiment. 

In  Section  6,  we  illustrate  the  application  of  the  statistical 
techniques  to  an  actual  queueing  simulation.  For  Illustration 
purposes,  we  choose  the  single-server  queue  with  exponential  inter¬ 
arrival  and  service  time  distributions,  since  its  theoretical 
properties  are  well-known  and  provide  a  means  of  comparison  with 
simulation  results.  Confidence  intervals  are  obtained  for  various 
3teady-state  quantities  of  interest,  including  the  mean  and  standard 
deviation  of  the  waiting  time,  the  expected  number  of  customers  in 
the  system,  the  expected  value  of  a  penalty  function  on  the  waiting 
time,  the  tail  of  the  waiting  time  distribution,  the  expected  number 
of  customers  served  in  a  busy  cycle,  the  expected  value  of  a  penalty 
function  or.  the  number  of  customers  served  in  a  busy  cycle,  the 
tail  of  the  distribution  for  the  number  of  customers  served  in  a 
busy  cycle,  and  others.  The  two-step  procedure  is  then  illustrated 
for  the  mean  waiting  time.  In  a  single  realization  of  a  second- 
step  run  consisting  of  18,000  busy  cycles,  a  98%  confidence  interval 
[.08£,  .108]  is  obtained  for  a  queue  with  an  arrival  rate  of  5 
customer^  per  time  unit  and  a  service  rate  of  10  customers  per  time 
unit  (the  theoretical  mean  waiting  time  is  0.1).  Two  systems  with 
different  service  rates  are  then  compared  by  means  of  a  confidence 
interval  on  the  difference  of  the  mean  waiting  times.  Finally,  a 
method  is  illustrated  by  which  one  may  infer  sensitivity  to  an 
unknown  parameter.  In  particular,  if  the  arrival  rate  is  unknown 
and  one  observes  the  system  at  two  different  values  for  the  arrival 
rate,  then  under  certain  linearizing  assumptions,  one  may  make 
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confidence  interval  statements  for  any  arrival  rate  between  the  two 
observed  values.  Indeed,  one  may  obtain  a  confidence  band  about  the 
expected  waiting  time  (or  some  function  of  the  waiting  time) 
expressed  as  a  function  of  the  unknown  parameter.  It  is  felt  that 
such  a  capability  could  be  an  effective  aid  in  assessing  the 
adequacy  of  the  input  data  accuracy,  and  in  efficiently  uncovering 
those  parameters  that  warrant  more  in-depth  analysis.  Of  course 
this  technique,  together  with  all  of  the  techniques  illustrated  in 
Section  6,  is  only  feasible  because  of  the  busy  period  structure 
which  is  exploited. 

Finally,  in  Section  7,  we  compare  our  approach  for  estimating 
the  mean  steady-state  waiting  time  with  a  commonly  used  alternative. 

In  our  approach,  one  arrives  at  an  estimate  essentially  by  forming 
the  sample  mean  of  the  customer  waiting  times  in  some  fixed  number 
of  busy  cycles.  The  actual  number  of  customers  observed  is  thus 
random.  In  an  alternative  approach,  one  might  simulate  a  fixed 
number  of  customers,  disregard  an  initial  fraction  of  the  customers 
as  transient  observations,  then  form  a  sample  mean  with  the  remaining 
customers.  In  both  cases,  the  estimate  obtained  converges  to  the 
theoretical  value  as  the  number  of  observations  becomes  large.  In 
our  approach,  however,  one  is  able  to  make  rigorous  statistical 
statements  by  taking  advantage  of  the  independence  of  the  busy  cycles. 
In  Section  7,  we  seek  to  determine  whether  this  advantage  is  obtained 
at  the  expense  of  reduced  accuracy  in  the  estimator.  Experimental 
results  are  obtained  for  the  single-server  queue  used  in  the  previous 
illustrations.  It  is  found  that  for  this  queue, statistics  which  are 
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obtained  in  costparahl.  length  computer  run.  for  the.e  two  approach,, 
are  very  clo..  in  accuracy.  Thu,,  one  My  apparently  „ak,  u»  of 
the  methods  of  this  report  vithont  fear  of  reduced  accuracy. 

Continuing  reeearch  in  this  area  will  extend  and  illustrate 
this  approach  for  more  general  ai.ul.tion.  and  describe  In  detail 
thft  steps  necessary  for  implementation. 


A  NEW  ATP ROACH  TO 


SIMULATING  STAI3LE  STOCHASTIC  SYSTEMS.* 

I  -  GENERAL  MULTI-SERVER  QUEUES 

by 

Michael  A.  Crane  and  Donald  L.  Iglehart 

1.  INTRODUCTION 

The  principal  goal  of  most;  simulations  of  stable  stochastic 
systemr.  is  to  estimate  propertieo  of  the  stationary  or  steady-state 
behavior  of  the  system.  Two  of  the  major  problems  in  such  simulations 
are  rhe  statistical  dependence  between  successive  observations  and 
the  inability  of  the  simulator  to  begin  the  system  in  the  9teady-state . 
The  first  problem  has  necessitated  using  methods  cf  time  series 
analysis  rather  than  classical  statistics.  The  second  has  inspired 
many  simulators  to  let  the  system  run  for  a  sufficient  length  of 
time  so  that  the  initial  transient  wears  off  and  a  steady-state 
condition  obtains.  This  procedure,  of  course,  requires  a  judgement 
on  how  long  to  let  the  system  run  before  making  observations. 

For  many  stochastic  systems  being  :,imulated  it  is  possible  to 
find  a  random  grouping  of  observations  which  produces  independent 
identically  distributed  (i.i.d.)  blocks  from  the  start  of  the 
simulation.  This  grouping  then  enables  the  simulator  to  avoid  the 
two  problems  mentioned  above.  He  has  at  uis  disposal  the  methods  of 
classical  statistical  analysis  such  as  confidence  intervals, 
hypothesis  testing,  regression,  and  sequential  estimation  since  the 
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observations  are  now  i.i.d.  Furthermore,  information  that  is  useful 
in  estimating  the  steady-state  behavior  of  the  system  can  be  collected 
from  rcratch  thus  eliminating  the  pioblem  of  the  initial  transient. 

The  key  requirement  for  obtaining  these  i.i.d.  blocks  is  that 
the  system  being  simulated  return  to  a  single  state  Infinitely 
often  and  that  the  mean  time  between  such  returns  is  finite.  This 
requirement  will  be  met  fer  many,  but  not  all,  stable  systems  that 
might  be  simulated. 

In  this  paper  we  shall  illustrate  the  main  ideas  of  this 
approach  in  the  context  of  the  GI/G/s  queue,  where  a  >_  1.  Other 
stochastic  systems  which  can  be  dealt  with  in  this  same  manner  will 
be  discussed  in  future  publications.  This  paper  is  organized  as 
follows.  Section  2  reviews  some  revelant  results  from  the  theory  of 
statistical  Inference.  Section  3  summarizes  the  probabilistic 
structure  of  the  GI/G/1  queue  with  an  eye  toward  using  these  results 
in  carrying  out  a  simulation.  In  Section  4  a  similar  treatment  is 
given  to  the  GI/G/s  queue  for  s  >  1.  In  Section  5  we  discuss  a 
two-step  procedure  for  estimating  the  stationary  expected  value  of 
a  general  function  of  the  waiting  time.  Numerical  illustrations 
for  the  M/M/1  queue  are  given  in  Section  6.  Finally,  in  Section  7 
we  co.xpare  our  procedure  with  a  commonly  used  alternative. 

2.  PRELIMINARY  RESULTS 

Before  proceeding  further,  it  is  useful  to  review  some  relevant 
ideas  from  the  theory  of  statistical  inference.  As  usual,  we  have  a 
given  probability  triple  and  a  sequence  of  random  variables  (r.v.'s) 
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defined  on  It.  Since  Che  detailed  construction  of  this  triple  is 
not  important,  we  omit  it  from  the  discussion.  Let  6  be  some 
unknown  parameter  and  let  £  tnd  6  be  two  random  variables 
obtained  from  a  statistical  experiment.  If  0  <Y  <  1  and 

P{£  £  0  <_  6}  £  1  -  y,  we  say  that  [£,  5]  is  a  confidence  interval 

for  0  with  at  least  1Q0(1-Y)X  confidence.  Roughly  speaking,  if 
many  values  are  obtained  for  £  and  3  in  independent  replications 
of  the  experiment,  the  interval  [£,  5]  will  surround  0  at  least 
100(1-y)X  of  the  time. 

Given  confidence  Intervals  for  one  or  more  parameters.  It  Is 
possible  to  obtain  confidence  intervals  for  various  functions  of 
these  parameters.  For  example,  suppose  0j  and  ©2  are  two 

unknown  parameters,  and  (£^,  5^]  and  f£2,  5^]  are  confidence 

intervals  for  0^  and  02  with  at  least  100(1-Y^)X  and  100(1-Y2)2 
confidence,  that  ia  P(£^  £  9^  £  8^}  —  ^_Y1*  and 
P{£2  £  ©2  £  ~  ^en  by  a  straightforward  argument, 

P{£t  "  ®2  —  91  ~  62  -  ®1  ~  -2^  -  l-Yi~Y2 •  (1) 

If  A  >  0,  then 

P{A£X  +  B  £  A9X  +  B  £  A0X  +  B}  £  1  -  Yj_-  (2) 

If  9^  £  0  and  a  >  0,  then 

P((t£1J+)3  £  ej  £  ([e1J+)a)  £  l  -  y:.  (3) 
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where  (x]  ■  max  (Q,x).  If  Oj  >  0  and  0^  >  0,  then 


till  e  e 

p{_i —  <  -L  <  - _}  >  1  -  y  -  y. 


8, 


*2 


+  -  ‘  '1  '2' 


(*> 


Finally,  suppose  0 (X )  is  an  unknown  parameter  which  is  a  linear 
function  of  the  parameter  X.  Suppose  X^  <  X^  and 
P{8  <_  61>  >_  1  -  y  and  P(02  <_  6(X  )  <_  >  1  -  yr 

Then 


(X-X  )(0  -  6  ) 

p{e  + - A — £- — L.  <  e(x)  _c  e,+ - ± — 

l  a2  -  *i  1  A2  ~  A1 


(X-X1)(52  -  Qj) 


for  a!3.  X^  <_  X  <_  X2}  ^  1  -  -  Y2» 


(5) 


This  fact  is  potentially  useful  for  studying  the  sensitivity  of  a 

nonlinear  function  0(X)  to  X  over  some  small  interval  [X^,  X2] . 

Now  let  {XQ:n  ^  1}  be  a  sequence  of  i.i.d.  r.v.'s  with 

2  2 

finite  second  moments,  and  define  9  »  E{X^}  and  o  »  o  (X^l.  Let 
$  denote  the  distribution  function  for  a  normal  random  variable  with 
zero  mean  and  variance  one;  i.e.. 


$ (x)  -  j 


$>(OdE  -»  <  x  <  », 


where 


♦  U) 


1  -r/2 

7Z7e 


**oo  <  £  <  oa 


For  0  <  x  <  1,  define  zx  =  i  *(x).  Also,  define  the  sample  mean 
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and  sample  standard  deviation 


X<„>  s  i 


n 

j-i ] 


n»l,2, . . . , 


»(n)  =  Ir^r  t  <x<  '  X(n))2] 
j-1  2 


2 ,1/2 


i  n— 2 13  y  • « #  • 


Now  we  know  by  the  central  limit  theorem  that  the 


lio  P{/n  (X(n)  -  6)/o  £  z}  •  $(z),  -®  <  z  <  <*>. 
n-H» 


Furthermore,  It  may  be  shown  using  the  strong  law  of  large  numbers 
that  the 


P{lio  s(n)  *  o}  ■  1. 

n-x» 


It  then  follows  by  a  well-known  result,  cf.  CHUNG  (1968),  Theorem 
4.4.8,  that  the 


lim  P{/n(X(n)  -  9)/s(n)  £  z)  -  $(t),  -»  <  z  <  00 

n-«o 


Hence,  for  0  <  y  <  1  and  large  n, 


?{~zl-yj2  -  '  ©>/s<n>  <  *  1-Y 


or 


P(X(n)  -  zi_y/2® (n)/*^  £  ©  £  X(n)  +  z^^bM  !</n)  *  1-y, 
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giving  an  approximate  100(l-y)X  confidence  interval  for  0.  This 
method  will  be  used  for  analyzing  sequences  of  i.i.d.  r.v.'s 
introduced  in  later  sections. 

3.  PROBABILISTIC  STRUCTURE  OF  THE  GI/G/1  QUEUE 

Consider  now  a  GI/G/1  queueing  system  in  which  the  Oth  customer 

arrives  at  time  ■  0,  finds  a  free  server,  and  experiences  a 

service  time  v0.  The  nth  customer  arrives  at  time  t  and 
0  n 

experiences  a  service  time  v^ .  Let  the  Interarrival  times 

t  -t  ,  -  u  ,  n>l.  Assume  that  the  two  sequences  {v  :n  >  0} 
n  n-i  n  —  n  - 

and  {u  :n  >  1}  each  consist  of  i.i.d.  r.v.'s  and  are  themselves 
n 

independent.  Let  E{u  }  ■  >.  E{v  }  ■  w  \  and  p  ■  X/u  where 

n  n 

0  <  X,  u  <  Thus  u  (X)  has  the  interpretation  cf  the  mean 
service  (arrival)  rate.  The  parameter  p  is  called  the  traffic 
intensity  and  is  the  natural  measure  of  congestion  for  this  system. 

We  shall  assume  that  p  <  1,  a  necessary  and  sufficient  condition 
for  the  system  to  be  stable. 

The  principal  system  characteristics  of  interest  are  Q(t) , 
the  number  of  customers  in  the  system  at  time  t;  W^,  the  waiting 
time  (time  from  arrival  to  commencement  of  service)  of  the  nth 
customer;  W(t),  the  v;ork  load  facing  the  server  at  time  t;  B(t), 
the  amount  of  time  in  the  Interval  [0,t]  that  the  server  is  busy; 
and  D(t),  the  total  number  of  customers  who  have  been  served  and 
have  departed  from  the  system  in  [0,t]. 

Here  we  shall  review  the  basic  structure  of  the  GI/G/1  queue 
relevant  to  our  simulation  study.  For  a  comprehensive  treatment  of 
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these  and  other  results  for  the  GI/G/1  queue  see  IGLEHART  (1971).  To 


begin  the  analysis  of  the  process  {W  :n  >  0}  let  X  *  v  ,  -  u 

n  n  n-1  n 

and  set  Sq  "0,  SQ  ■  +. . .+  Xn»  n  >_  1.  The  following  recursive 

relationship  exists  for  the  Wn's: 


By  Induction  one  can  show  that 


W  ■  max{S  -S,  :k  ■  0,1,. ,.,n),  n  >  0. 
n  n  k  — 


Using  the  notion  of  optional  r.v.'s,  it  can  be  shown  that  there  exists 


a  sequence  of  r.v.'s  1.  0}  such  that  Bq  d  °’  Bk  *  9k+l  ,  and 


■  0  with  probability  one.  In  other  words,  the  customers  numbered 
are  those  lucky  fellows  who  arrive  to  find  a  free  server  and 


experience  of  no  waiting  in  the  queue.  The  fact  that  there  exists  an 


infinite  number  of  such  customers  is  a  direct  consequence  of  the 
assumption  that  p  <  1.  The  time  axis  •  t0,“>)  can  be  divided 
into  alternating  intervals  during  which  the  server  is  busy,  idle,  busy, 
etc.  We  call  these  intervals  busy  periods  (b.p.'s)  and  idle  periods 
(l.p.'s).  An  i.p.  plus  the  preceding  b.p.  is  called  a  busy  cycle 


(b.c.).  If  we  let  •  6^  -  k  21  1>  then  represents  the 

number  of  customers  served  in  the  kth  busy  period  (b.p.)  and  they 
are  numbered  ^-1  +  1, . . .  ,6j.-l)  • 
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Next  define  the  random  ve  :tora 


{V  -6.  1+1,,,,’\}  *  k  -  lm 

k-1  k 


*k"(VrV  and 

Note  that  these  vectors  are 


unusual  In  the  sense  that  they  have  a  random  number  of  components , 
namely,  a^+1.  Observe  that  the  vector  -  {a^,  X^,...,Xa  } 
Includes  all  the  data  required  to  completely  construct  the 
behavior  of  the  system  in  the  first  b.p.  The  principal  fact  that 
permits  us  to  decompose  the  system  into  1.1. d.  blocks  is 
that  the  Y^'s  are  i.i.d.  Hence  we  have  the  intuitively  plausible 
conclusion  that  comparable  r.v.'s  in  different  b.p.'s  are  i.i.d. 

Now  define  the  following  r.v.'s  for  k  ^  1: 


v  /■••+  x-i’ 

k-1  k 


ck  "  uek_1+iK'-'KV 


\  “ck  "  V 


YkX)  -  i  (Sg  +j 

K  J-0  Ck-1  3 


-  Sc  ), 


6, 


k-1 


*k  xk  nk’ 


ana 


V1 


r»>  -  r  ((s 


j-0 


Sk-1+J '  Vi*  Vsk-l+3  +  *  \-i+sh 


These  r.v.'s  have  the  following  interpretations:  is  the  length 

of  the  kth  b.p.,  £k  the  length  of  the  kth  b.c.,  the 

length  of  the  kth  i.p.,  Y^^  the  sum  of  the  waiting  times  in 
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(2) 

kth  b.p.,  Yk  the  integral  under  the  curve  Q(s)  in  the  kth 

(3) 

b.p.,  and  Y£  the  integral  under  the  curve  W(a)  in  the  kth  b.p. 

Because  the  V^'a  are  l.i.d.  so  are  the  r^’s,  5k's,  vit  8’  Yk^^'8’ 

C2>  (3') 

Yk  7,s,  and  Y^  7 'b.  It  is  these  i.i.d.  sequences  that  we  shall 

observe  in  our  simulation  with  an  eye  toward  using  classical 

statistical  methods  of  analysis. 

Next  ve  record  some  known  results  on  expected  values  of  these 

2 

r.v.'s.  We  shall  assume  that  the  E(Vq)  <  "  •  Then  the 

on 

Ela^}  -  exp{  £  n_1  P[Sn  >0]},  E{ "  u'1  E{ak>, 

-i  -i  m 

E{tk>  -  X  X  E{ok),  E{vk>  «  X  A  (1-p)  EC^}  ,  E{Y^AM  -  E{W}  E^}  <  », 

E{Y*2))  -  (XE{W)  +  p)  E{£k),  and  E{y£3)}  -  (pE{W)  +  ±  XE{v2))  EU^}, 

where  W  is  the  r.v.  to  which.  Wq  converges  in  distribution.  In  other 

words,  W  is  the  so-called  stationary  waiting  time.  An  expression  for 

itS+ 

the  characteristic  function  of  W  exists,  in  terms  of  the  E{e  n}, 
however,  it  is  difficult  to  evaluate  for  specific  cases.  Also 

OD 

E(W)  “  H  n  1  E(Sn)  •  If  u^  has  a  non-lattice  distribution,  it 
n»l 

is  known  that  Q(t)  Q  and  W(t)  *^W*  as  t  ■+•  °°,  where  — >  denotes 
weak  convergence  (convergence  in  distribution).  Hence  Q  is  the 


stationary  queue  length  and  W*  is  the  stationary  virtual  waiting 


time.  Furthermore,  E { Q >  •  XE{W}-i-p  and  E(W*}  ■ 


which  are  the  terms  appearing  in  the  expressions  for 

E{Y^3)}  . 

k 


pE(W)  +  XE(v2>/2 

E{Y,(2)}  and 
k 


More  generally  let  f  be  a  measurable  function  from  [0,®)  i  ito 


(-oo, os).  Assuming  that  the  E{f(W)}  <  ®  ,  then 

V1 

E{  £  f (Sg  +j  “  Sg  )  }  -  E{ f (W) )  E{ak)  •  This  identity  plus  the 

j-0  k-i 

fact  that  the  Vk's  are  allows  us  to  treat  E(f(W))  in  the 
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game  manner  that  we  do  E{W}.  In  Section  6  a  variety  of  functions  f 
of  practical  Interest  will  be  discussed. 

While  the  parameters  E{a^},  P.{W},  E{Q),  E{W*},  and  E{f(W)} 

can  theoretically  all  be  calculated  from  the  distributions  of  Vq  and 
u^,  these  calculations  are  very  difficult  and  one  might  choose  to 
estimate  them  through  simulation.  It  is  problems  of  this  sort  that  we 
address  in  this  paper. 

4.  PROBABILISTIC  STRWURE  OF  THE  Cl/G/a  QUEUE 

Now  suppose  there  are  s  >  1  servers.  The  appropriate  value 
of  the  traffic  intensity  is  now  p  »  A/oU  ,  which  ve  assume  again  is 
strictly  less  than  one.  It  is  known  (see  KIEFER  and  WOLFOWITZ  (1955, 
1956),  and  LOYNES  (1962))  that  p  <  1  is  a  necessary  and  sufficient 
condition  for  the  queueing  system  to  be  stable.  However,  p  <  1 
is  not  enough  to  Insure  that  the  queue  becomes  idle  infinitely  often, 
which  is  the  key  to  the  analysis  in  the  case  s  *  1.  A  simple  and 
reasonable  condition,  which  we  shall  henceforth  assume,  was  recently 
provided  by  WHITT  (1971);  namely,  that  the  *'(un  >  v^_^}  >  0* 

Additional  results  on  this  problem  are  contained  in  KENNEDY  (1971). 

The  analysis  of  this  queue  is  best  approached  through  a  vector¬ 
valued  waiting  time  process  {W  :n  >,0}  introduced  by  KIEFER  and 
WOLFOWITZ  (1955,  1956).  We  assume  the  customers  are  served  by  the 
first  available  server.  If  more  than  one  server  is  free,  the 
customer  is  served  by  that  server  with  the  smallest  index.  For  a 
fixed  realization  of  the  queueing  system,  think  of  assigning  the 
customers  when  they  arrive  to  the  server  who  will  eventually  serve 
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them.  Let  W  •  (W  W  )  be  the  vector-valued  process  whose 

-n  nl  ns  r 

components  W  ^  denote  the  workload  of  the  server  with  the  1th 
lightest  load  just  prior  to  the  arrival  of  the  nth  customer.  Thus 
is  the  actual  waiting  time  of  the  nth  customer.  The 
sequence  can  be  generated  recursively  as  follows: 


W  ,,  -  +  y  ,  -  U  )]+. 

'•n+l  1  'Dn  n-1  *n  1  J 

where  V  ,  -  (v  , ,  0, .  . . ,0) ,  U  -  (u  ,...,u  ),  [X]+  -  <xt,...,x*) 
-n-i  n-i  -n  n  n  i  s 

S  SB 

for  X  ■  (x^,...,xa)  e  (R  ,  and  F:&  -*•  (R  rearranges  the  components 

of  its  argument  in  ascending  order.  It  is  clear  from  this 
representation  of  Wr  that  {\Jn:n  1  0)  is  a  Markov  process  with 
state  space  E  ■  (x  e  R8 :0  ±  x^,  •  •  •  £  xg).  The  condition  of  Whitt 

mentioned  above  guarantees  that  the 

P{Wn  ■  0  infinitely  often}  *  1 

and  that  the  expected  time  between  visits  to  0  is  finite.  This 
condition  allows  us  to  define  the  sequences  of  r.v.'a  (a^:k  >  1}  , 
{B^ik  0},  and  {V^tk  >_  1}  just  as  in  the  case  of  the  single- 
server  queue.  In  this  case  customers  numbered  6^.  arrive  to  find 
all  servers  idle. 

Thus  we  can  again  proceed  to  decompose  the  queueing  system 
into  Independent  b.c.'s  just  as  was  done  in  the  case  s  *  1.  Now 
a  b.c.  is  defined  as  an  interval  of  time  during  which  at  least  one 
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server  is  busy.  In  fact  we  can  take  over  the  same  notation.  Hoveypr, 


in  this  case  we  have  no  closed  expression  for  E{c<k}  and  thus  we  have 
no  alternative  but  to  carry  out  a  simulation.  Again  F, { )  “  h  ^  ECc^}, 


EUk) 


X-1E{uk),  E{u  }  -  A_1(l-p)  E{ak),  and  EfY^15}-  E{W)  Efc^.} 

r(l) 


where  W  is  again  the  stationary  waiting  time  and  Yk  Is  the  sum  of 
the  waiting  times  in  the  kth  b.p.  There  is  no  closed  expression  for 
ElW}  and  again  simulation  seems  like  the  only  recourse. 


5.  A  TWO-STEP  PROCEDURE  FOR  ESTIMATING  Eff(W)} 

Let  f  be  some  positive  measurable  function  on  the  customer 
waiting  time  and  let 


Y 


k 


k-1 


). 


We  know  from  the  general  theory  that  {Yk:k  >  D  is  a  sequence  of 
l.i.d.  r.v.'s  and  that  EfY^  -  E{f(W)}  Efc^}.  This  last  relation 

is  at  the  heart  of  the  method  to  be  proposed.  By  observing  the 
sequences  {Yk;k  1}  and  {c^tk  _>  1}  and  obtaining  a  100(l-Yj)% 
confidence  interval  for  E{Y^}  and  a  100(1-Y2)^  confidence  interval 
for  E{a^},  it  is  possible  to  obtain  a  100(1-y^  -  Y2^  confidence 
interval  for  E{f(W)}  using  (4).  Two  questions  naturally  arise: 
given  a  desired  100(1-y)%  confidence  interval  for  E{f(W)},  how 
should  y^  and  Y2  chosen,  and  what  effect  do  these  choices 
have  on  the  required  sample  size  N  ?  It  is  these  questions  which 
we  address  in  this  section.  Given  a  desired  y»  we  develop  a  two- 
step  procedure  in  which  the  simulator  chooses  y^»  Y2»  and  N 

after  obtaining  preliminary  estimates  for  system  parameters  in  an 
initial  "short"  simulation  run. 


L2 


The  full  procedure  consists  of  two  independent  simulation  runs, 
the  first  for  a  total  of  m  busy  cycles  and  the  second  for  N  busy 
cycles,  where  N  is  to  be  determined.  The  initial  run  is  a  short 
run  which  serves  an  a  planning  guide  for  the  much  longer  second  run. 
The  first  run  provides  the  simulator  with  estimates  of  confidence 
interval  lengths  for  E{f(W)}  which  could  be  obtained  in  the  second 
run  for  various  choices  of  V2>  am*  N.  For  each  N,  the 

estimated  length  is  minimized  over  all  possible  values  of  and  Yj 

with  +  y2  “  ant*  m:1-n^na^  length  is  expressed  as  a  function 
of  N.  The  simulator  is  thus  able  to  consider  tradeoffs  between 
the  length  of  confidence  interval,  level  of  confidence,  and  cost  of 
computer  resources. 

Let  (Y,,...,Y  )  and  (a, ,...,a  )  be  observations  during  the 
i  m  J-  tn 

initial  run,  and  define 


Given  these  statistics,  the  simulator  can  estimate  the  size  of 
confidence  intervals  for  E{Y^}  and  E{a^}  which  could  be  obtained 
in  a  second  run,  for  various  values  of  N,  y^,  and  v0.  If  the 
second  run  were  to  be  made  with  sample  size  N,  -  x,  and  Yj  *Y  - 
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then  estimates  for  lOOfl-y^JZ  and  100(1-y2>%  confidence  intervals  for 
E{Y^}  and  E{a^}  are 


[Y  -  z l  x/2  yrf.  Y  +  z1-x/2  Sy/*^ 


and 


[a  -  Z 


l-y/2  +  x/2 


sa/^. 


a  +  z 


l-y/2  +  x/2 


so//N], 


Using  (4),  an  estimate  for  the  length  of  a  100(1-y)X  confidence 
interval  for  E{ f (W) }  is  therefore 


.  ll  ~  _ 

l“-*l-,/2+,/28«/'5l+  “+2l-Y/2+»/28a/^ 


For  a  fixed  N  and  y,  it  is  desirable  to  choose  that  x  which 
minimizes  H(N,y,x).  Using  the  expansion 


A-i  E  <-»3  fc1  • 


3-0 


we  have,  for  large  N,  the  approximation 


2s. 


2Y  s 


i(N.Y.x)  *  F(N ,y*x)  -  —  Zl_x/2  +  Z^T  *!.y/2  +  x/2-  <6> 


We  thus  consider  the  problem  of  choosing  x  so  as  to  minimize 
F(N,y,x),  subject  to  0  <  x  <  y- 
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Recalling  that  m  i  ^(y)  for  0  <  y  <  1,  it  may  be  easily 
shown  that  F(N,y,x)  is  a  strictly  convex  continuous  function  of  x 
for  0  <  x  <  y  and  that  F(N,y,x)  *  as  x  -►  0  or  x  -*•  y.  It 

follows  chat  a  unique  minimum  is  obtained  at  x  -  x^,  the  unique 

£ 

root  in  (0,y)  of  the  equation  —  F(N,y,x)  ■  0. 

oX 

Now 


3^  F(N*y*x) 


_ 1 _ 

*<r1a-x/2)) 


i _ 

<K$”1(l-y/2+x/2)) 


</2va^  j  2  >/^1r^8a  1  2 

“  Hf  CXP{2  Z1-y/2+x/2)’ 


Setting  this  equal  to  zero,  we  obtain 

2  2  a®Y 

Zl-Y/2*x0/2  "  Z1-xq/2  "  Un^~  2  C>  (?) 

a 

Note  that  the  solution  Xq  obtained  in  (7)  is  independent  of  N. 
That  is,  for  a  fixed  y,  there  is  a  unique  choice  y^  -  x^  which 
ia  optimal  for  every  value  of  N,  The  estimate  of  the  length  of  the 
confidence  interval  for  Elf(W)}  is  obtained  from  (6): 


where 


£q(N,y)  -  i(N,y,x0)  a  D(y)/^?5 


(8) 


D(y) 


2s„  2Ys 

_ Y  +  _ a 

S  *W2  =  Zl-Y/2  +  xft/2- 


(9) 


Given  (8)  end  (9),  the  simulator  can  make  a  rational  choice  of  N 
and  y  based  on  a  consideration  of  computer  time  available,  desired 
level  of  confidence,  and  desired  length  of  confidence  interval. 

Table  1  shows,  for  several  values  of  y,  the  quantity 
2  2 

zl_y/2  +  x/2  ~  zl-x/2  an  a  f,-,nctio‘'1  °*  Y  and  h  =  x/y.  For  a 
given  value  of  c,  computed  after  the  initial  simulation  run,  and  of 
Y  ,  Table  1  may  thus  be  used  to  find  that  value  x^  solving  (7). 

D(y)  can  then  be  computed  from  (9)  using  standard  tables  for  the 
normal  distribution.  These  computations  are  illustrated  in  Section  6. 

It  should  be  noted  that  (8)  and  (9)  give  only  an  estimate  of 
the  length  of  the  confidence  interval  which  could  be  obtained  in  th»* 
second  simulation  experiment  of  N  busy  cycles,  this  estimate  being 
made  at  the  end  of  the  first  experiment.  The  actual  confidence 
interval  is  computed  at  the  end  of  the  second  experiment,  using 
confidence  intervals  for  E(Y^}  and  E{a^}  together  with  (4). 

It  is  interesting  to  contrast  the  method  outlined  above  with 
classical  sequential  estimation  procedures*,  cf,  ANSCOMBE  (1953). 

In  the  classical  procedures,  one  generally  desires  to  compute  a 
confidence  interval  of  a  fixed  length  for  a  single  unknown  parameter. 
The  interval  length  and  confidence  level  are  fixed  at  the  start  of 
the  procedure,  and  the  number  of  observations  required  is  then 
dictated  by  the  experimental  observations  and  a  sequential  stopping 
rule,  beyond  the  control  of  the  experimentor .  In  contrast,  the 
two-step  procedure  outlined  in  this  section  deals  with  estimation 
of  the  ratio  of  two  unknown  parameters,  so  that  one  must  consider 
the  problem  of  minimizing  the  interval  length  for  the  ratio  over  all 
possible  lengths  for  the  individual  parameters.  Furthermore,  any 
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decisions  involving  the  desired  degree  of  confidence,  the  Interval 
length,  and  the  ultimate  number  of  observations  are  deferred  until 
the  end  of  the  first-step  experiment,  when  preliminary  estimates  for 
the  parameters  are  available.  The  first  step  thus  provides  a  basis 
for  the  economic  decisions  which  must  be  made  by  the  simulator  in 
using  limited  computer  resources. 


6.  SOME  NUMERICAL  ILLUSTRATIONS  FOR  THE  M/M/1  QUEUE 

In  this  section,  we  give  numerical  examples  which  combine  the 
statistical  techniques  of  Sections  2  and  5  with  the  queueing  results 
of  Sections  3  and  4.  The  M/M/1  queue  is  used  for  illustration 
because  its  theoretical  properties  are  well-known  and  provide  a  means 
of  comparison  with  the  simulation  results.  In  this  case,  the  inter¬ 
arrival  and  service  times  are  distributed  as  exponential  r.v.'s  with 
parameters  A  and  y  respectively. 

From  Section  3,  we  know  that  the  sequences  {cc^sk  >  1 } , 

{nk:k>l},  Uk:k>l},  (vk:k  >  1),  (Y^k  >  1],  {Y^2):k>l}, 

(3) 

and  {Yk  :k  1)  each  consist  of  i.i.d.  r.v.'s  and  are  thus 
amenable  to  statistical  analysis.  In  order  to  demonstrate  the  power 
and  versatility  of  our  methods,  we  define  five  additional  sequences 
of  i.i.d.  r.v.'s.  For  k  >  1  let 


“k"1 

'(4)  -  V  f  (W 
K  j-0  4 


> 5  i  (v^2' 


,<5> . 


v1 


V1 

I  £5(W6  +j>  E  I  [(w 

j-0  3  Bk-1  J  j-0 


\-l+J 


•1)+]1/2, 
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i 


\ 


k 


\ 


l 


,<6) 


V1  V1 

^6^8  +1^  =  ^{W  >  2}  * 

j-0  0  ^k-l  3  j-0  lw6k_j+J  -  *  1 


rO)  m 


o,  a  -  1,2 

f.(o  )=<  5,  aj  -  3,4,5 

10,  otherwise, 


,(8)  . 


f8(ak) 


«J  lj  “k-1*2 


0, 


otherwise, 


where  the  indicator  function  ^  takes  on  the  value  1  on  the  set  A 

and  0  otherwise.  The  function  f^  Is  usef’  l  for  estimating  the 

second  moment  and  variance  of  tie  stationary  wai  h.fi  time  W,  since 
2 

E(f4(W)}  ■  E{W  }.  The  function  f^  can  be  interpreted  as  a 

penalty  function  on  the  waiting  time  of  the  customer.  The  function 
fg  provides  a  means  for  estimating  the  tail  of  the  stationary 
waiting  time  distribution,  since  E(fg(W)}  ■  P(W  .2).  The 
function  ty  can  be  interpreted  as  a  penalty  or  cost  function  on 
the  number  of  customers  served  during  a  busy  period.  This  might  be 
appropriate,  for  example,  if  servers  are  to  be  relieved  at  the  end 
of  each  busy  period  and  a  penalty  must  be  paid  when  a  server  is 
required  to  serve  an  excess  number  of  customers  without  relief. 
Finally,  the  function  fg  provides  information  on  the  distribution 
of  a^,  since  E{fg(a^)}  -  P{a^  -1  or  ■  2}. 

In  order  to  illustrate  the  computation  of  confidence  intervals 
for  the  various  parameters  of  interest,  a  single-step  simulation 
run  was  implemented  with  X  ■  5,  p  ■  10,  and  N  -  2000  busy 
cycles  (demonstration  of  the  two-step  procedure  of  Section  5  follows 
later).  That  is,  2000  observations  were  made  from  each  of  the 
sequences  given  above. 


i 
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Now  we  have  seen  earlier  that  if  is  a  sample  of 

i.i.d.  r.v.'s  with  sample  mean  X,  sample  standard  deviation  et  and 
finite  second  moments,  then  [X  -  zy_yj2  X  +  zi-y/2  *8  an 

approximate  100(l-y)2  confidence  interval  for  E{X^}  for  large  N. 

In  this  manner,  we  obtain  approximate  confidence  intervals,  shown  in 

Table  2,  for  Efc^},  Efr^},  E{^},  Etv^,  and  EfY^n)},  n  -  1,2 . 8 

based  on  the  observed  samples  of  2000  (The  moment  condition  is  satisfied 
because  all  of  the  moments  of  and  are  finite.)  Given 

approximate  confidence  intervals  for  E(Y^},  E{Y^},  E{Y^},  E{Y^}, 


E{a, },  we  use  (A)  to  obtain 

approximate  confidence  intervals  for 

E{W) ,  E(f4(W)),  E{ f j (W)  } , 

and  E{f, (W) } ,  since 

0 

eCyJ^} 

e(w} 

^  1 

ECa^} 

e(y1C4)  } 

Efa^ 

E(f4(W)} 

E(f5(W)} 

e{y<:>)} 

E(a1 } 

and 

E{Y(6)) 

E(f6(W)} 

1 

E{a.  ) 

Given  an  approximate  confidence  interval  for  E{W},  we  use  (2)  to 
obtain  an  approximate  confidence  interval  for  E{Q)  *  AE{W)  +  p. 

Finally,  given  approximate  confidence  intervals  for  E{W}  and  for 
2 

E(f4(W)>  ■  E(W  },  we  use  (1)  and  (3)  to  obtain  an  approximate 

2  2  1/2 

confidence  Interval  for  the  standard  deviation  a{W}  ■  [E{W  }  -  (E{W})  J  '  . 


and 
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Table  2  summarizes  the  results  of  the  simulation  run.  For  each 
(8) 

entry  through  E{Y^  '},  the  value  given  in  the  table  for  the  "point 
estimate"  is  the  value  of  the  sample  mean.  For  each  of  the  remaining 
entries,  the  point  estimate  is  the  appropriate  function  of  the  other 
point  estimates;  e.g.,  the  estimate  for  E{W}  is  the  ratio  of  the 
estimates  for  E{Y^}  and  E{a^}. 

Table  3  shows  the  point  estimates  and  90X  confidence  intervals  for 
E{W}  in  ten  replications  of  the  experiment.  As  can  be  seen,  for  these 
runs,  each  of  the  ten  intervals  surrounds  the  true  mean  E(U)  *  .1, 

We  next  demonstrate  the  two-step  procedure  of  Section  5  for 
obtaining  a  confidence  interval  for  E(f(W)}  .  For  illustrative 
purposes,  we  let  f  be  the  identity  function,  so  that  we  estimate 
E{W} . 

In  order  to  illustrate  the  procedure  a  first-step  run  of  500 
busy  cycles  was  made.  Estimates  obtained  from  this  run  are 

Y  -  .178  , 

a  -  1.95  , 

Sy  -  .749  , 

s  -  2.20  , 

a 

so  that 

asY 

<;  ■  2  In  *  2.63. 

a 

The  optimal  value  y^  =  x^  may  now  be  obtained  from  Table  1.  For 
several  values  of  y.  Table  4  shows  the  optimal  y^,  y^  ■  y-y^, 
and  the  associated  D(y).  For  a  second-step  simulation  of  N  busy 
cycles,  Iq(N,y)  “  D< y) / is  the  estimate  for  the  minimum  length 
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TABU  2 


SIMULA T ION  RESULT.''  FOR  THE  M/M/1  QUEUE 

(arrival  rate  X  ■  5 ;  service  rate  u  ■  10;  N  -  2000 
observed  busy  cycles) 


Parameter 

Theoretical 

Point 

Confidence 

Level  of 

Value 

Estimate 

Interval 

Confidence 

Et^) 

2.000 

2.110 

[1.994,  2.226] 

95% 

E  { > 

0.200 

0.215 

[.199,  .231] 

95% 

E^} 

0.400 

0.416 

(.396,  .435] 

95% 

E(v1) 

V1 

0.200 

0.201 

{.192,  .210] 

95% 

e{y[1} )  - 

tit  V 

j-o  -1 

0.200 

0.232 

[.187,  .277] 

95% 

e(yJ2))  - 

E  {  f  Q(s)da  } 

0.40Q 

i 

0.447 

[.387,  .507] 

95% 

E{Y^3)  )  - 

E {  /  W(s)ds ) 

jo  , 

i 

0.040 

1 

0.045 

[.038,  .053] 

95X 

e{y*4)  }  - 

V1 

E{£  (W  )2) 

J-0  3 

0.080 

0.096 

[.066,  .127] 

95% 

E{Y*5)  }  - 

“l"1 

e(£  /(W  -.1)  ) 

j-0  3 

0.240 

0.280  ( 

[.226,  .333] 

95% 

e(y£6)  >  - 

aI_l 

E%  I{Wj>.2) 

0.368 

0.438  ; 

(.357,  .519] 

95% 

e(y*7)  )  - 

E{f4(«1)  ) 

1.225 

1.375 

[1.248,  1.502] 

95% 

e(y<8>  )  - 

P(a^«l  or  a^«2  } 

0.815 

0.792 

[.774,  .810] 

95% 

E(W) 

0.100 

0.110 

[.084,  .139] 

90% 

E  { r 4 CW)}  " 

E(W2} 

0.040 

0.046 

[.030,  .064] 

90% 

E{f 5(W) )- 

E{/(W-.l)r) 

0.120 

0.133 

[.102,  .167]  | 

90% 

E{fg(W))- 

P(W>_.2) 

0.184 

0.208 

[.160,  .260]  | 

90S 

E(Q ) 

1.000 

1.050 

(.920,  1.195) 

90% 

o{W) 

0.173 

0.182 

(.101,  .238] 

80% 
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ESTIMATES  FO 
5*.  u  -  10;  N 


Replication 


TABLE  3 


[W}  IN  TEN  SIMULATION  REPLICATIONS 

)0  observed  busy  cycles;  level  of  confidence  >  90X) 


Point  Estimate 

Confidence  Interval 

0.110 

(.084,  .139] 

0.091 

1.071,  .114] 

0.095 

[.075,  .117] 

0.111 

(.075,  .151] 

0.096 

[.073,  .122] 

0.100 

[.077,  .126] 

0.092 

[.071,  .116] 

0.099 

(.074,  .128] 

0.096 

[.073,  .122] 

0,090 

[.068,  .115] 
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TABLE  4 


VALUES  FOR  y x .  Y2>  AND  D(y)  AS  A  FUNCTION  OF  y 
(Based  on  a  simulation  run  with  X  *  5,  u  “  10,  m  ■  500  observed  busy  cycles.) 


Y 

Y1 

Y2 

D(y) 

.2 

.167 

.033 

1.50 

.1 

.082 

.018 

1.83 

.05 

.041 

.009 

2.11 

.02 

.0164 

.0036 

2.44 

.01 

.0081 

.0019 

2.67 

i 


TABLE  5 

SIMULATION  RESULTS  FOR  THE  TWO-STEP  PROCEDURE 
(  X  •  5,  u  ■  10,  N=  18,000  observed  busy  cycles) 


Parameter 

Theoretical 

Value 

i 

Point 

Estimate 

Confidence 

Interval 

Level  of 
Confidence 

E{ax} 

2.000 

1.983 

[1.932,  2.035] 

99.64* 

E{Yj1} } 

0.200 

0.192 

[0.175,  0.209) 

98.36* 

E{W) 

0.100 

0.097 

[0.086,  0.108] 

98.00* 

confidence  interval  for  E(W) ,  and  Figure  1  shows  1q(N,y)  as  a 

function  of  N  and  y.  For  example  if  N  is  chosen  to  be  4000 

(about  a  5  second  run  on  the  IBM  360/67)  and  y  -  .05,  then  the 

estimated  confidence  Interval  length  is  .033,  for  a  deviation  of 

±18%  about  the  point  estimate  for  E(W}. 

Using  Figure  1  as  a  guide,  it  was  decided  that  the  second 

experiment  be  run  with  N  -  18,000  and  y  ■  .02  (98%  confidence). 

Table  5  shows  the  actual  results  of  this  run.  The  confidence  interval 

obtained  for  E(W)  is  of  length  0.022,  compared  to  a  length  of 

0.018  estimated  at  the  end  of  the  first  step. 

Suppose  now  chat  one  desires  to  make  comparative  Inferences 

about  two  different  queueing  systems .  For  example,  one  might  be 

interested  in  comparing  the  mean  stationary  waiting  times  E(W^} 

(2) 

and  E{WV  ')  for  systems  with  p  -  10  and  u  “  16,  each  having 

an  arrival  rate  X  -  5.  One  method  of  comparison  is  to  obtain  a 

confidence  Interval  on  the  difference  in  the  mean  waiting  times. 

This  can  be  done  by  computing  confidence  intervals  on  the  means 

E(W^}  and  E{W^}  separately,  then  using  (1)  for  the  difference. 

To  illustrate,  a  run  of  18,000  busy  cycles  was  made  with 

(2) 

u  “16,  resulting  in  a  98%  confidence  interval  for  E(W  )  of 
[.025,  .030]  (y^  and  Yj  were  chosen  the  same  as  earlier  for 

u  ■  10).  Combining  this  with  the  previously  computed  interval  for 
E(W^},  we  obtain  the  interval  [.056,  .083]  for  E{W^}  -  E{W^}. 
That  is ,  with  confidence  not  less  than  96%,  the  mean  stationary 
waiting  time  for  u  •  10  exceeds  that  for  p  ■  16  by  at  least 
.056  but  not  more  than  .083.  Statements  of  this  type  would  be  useful 
for  assessing  the  relative  value  of  different  proposed  system 


modifications. 
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Figure  1.  ESTIMATED  LENGTH  OF  CONFIDENCE  INTERVAL  FOR  l’HE  MEAN  WAITING  TIME  AS  A  FUNCTION  OF  NUMBER  0"  BUSY 
CYCLES  AND  LEVEL  OF  CONFIDENCE  (based  on  a  first-step  run  of  500  busy  cycles  with  arrival  rate 


Suppose  one  wishes  to  study  the  sensitivity  of  the  queueing 
system  over  a  range  of  parameter  values,  say  for  3  ^  ^  <^5  with 

V  -  10.  We  illustrate  again  for  the  mean  waiting  time  E{W(X)}. 

If  one  assumes  that  E{W(\)}  is  approximately  a  linear  function  of 
X  for  3  <_  X  <_  5,  then  (5)  may  be  used  to  obtain  an  approximate 
confidence  band  about  the  function  E { W( X ) }  over  that  interval. 

To  Illustrate,  a  run  of  18,000  busy  cycles  was  made  for  the  case 

X  -  3  and  u  *10,  and  a  98%  confidence  interval  obtained  for 

the  mean  waiting  time  is  (.037,  .045].  A.n  interval  [.086,  .108] 

was  previously  obtained  for  X  ■  5,  and  u  ■  10.  Then  under  the 

linear  assumption  for  E { W (X ) }  between  X  ™  3  and  X  »  5,  a 

96%  confidence  band  for  E{W(X)}  is  shown  in  Figure  2  together  with  the 

true  function  E{W(X)}.  That  is,  with  at  least  96%  confidence, 

.037  +  ( .049)  (X-3) /2  <_  E(W(X)}  _<  .045  +  (.063)  (X-3)/2  for  all 

3  <.  X  _<  5. 

The  same  method  could  be  combined  with  the  previous  method  for 
comparing  alternative  systems.  For  example,  if  E(W(X,p)}  is  the 
mean  waiting  time  for  a  system  with  parameters  X  and  p,  and 
if  vij  and  ^2  are  two  alternative  values  for  u,  one  might  wish 
a  confidence  band  for  g(X)  =  (E{W(X,u^)}  -  EfWCA,^)}]  over  a 
range  X^  _<  X  <_  X^.  Tha  is,  inferences  would  be  made  on  the 
sensitivity  of  the  system  differences  over  the  given  range  for  A  , 
based  ou  the  assumption  of  linearity  over  this  range.  This  is  done 
by  making  runs  at  only  four  different  parameter  settings:  (X^,  p^) , 

(X^vi  2 )’  (Xj.u^K  and  (Xj.l^)*  Figure  3  illustrates  our  results. 
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This  graph  illustrates  one  particular  experimental  realization,  based  on 
observations  at  the  two  endpoints  X  ■  3  and  X  =  5.  In  at  least  98%  of  such 
realizations,  the  confidence  limits  at  either  endpoint  will  surround  the  true 
value.  In  at  least  962  of  the  realizations,  the  confidence  limits  at  all 
points  3  <  X  <_  5  will  surround  the  true  values.  Thus,  for  example,  a 
98%  confidence  interval  at  X  -  4  Is  [.062,  .076]  in  this  realization. 


0  12  3  4  5 

Arrival  Rate,  X 


Figure  2.  A  962  CONFIDENCE  BAND  FOR  THE  MEAN  WAITING  TIME,  AS  A  FUNCTION  OF  THE 
ARRIVAL  RATE  A  OVF.R  3  <  X  ^5.  (Service  rate  p  -  10;  N  -  18,0 
observed  bu9V  cycles.) 


Difference  in  Mean  Waiting  Times,  E{W(X,u  )}  -  £{W(X, 


This  graph  Illustrates  one  particular  experimental  realization,  based  on 
observations  at  the  two  endpoints  X  ■  3  and  X  -  5.  In  at  least  96%  of 
such  realizations,  the  confidence  limits  at  either  endpoint  will  surround 
the  true  value.  In  at  least  92%  of  the  realizations,  the  confidence  limits 
at  all  points  3  <_  X  <5  will  surround  the  true  values.  Thus,  for 
example,  a  92%  confidence  interval  at  X  =  4  in  [.039,  .057)  in  this 


Figure  3.  A  92%  CONFIDENCE  BAND  FOR  THE  DIFFERENCE  IN  MEAN  WAITING  TIMES 
BETWEEN  TWO  SYSTEMS  WITH  DIFFERENT  SERVICE  RATES,  AS  A  FUNCTION 
OF  THE  ARRIVAL  RATE  X  OVER  3  <  X  <5.  (Service  rates 
y.  »  10,  ii.  ■  16;  N  =  18,000  observed  busy  cycles.) 


The  key  assumption  la  of  course  the  linear  relationship  over 
the  range  of  Interest.  This  ao3umptlon  Is  fairly  accurate  In  the 
c::?mples  given,  as  seen  It?  Figures  2  and  3.  In  future  work,  wc  Intend 
to  explore  the  possibility  of  obtaining  confidence  bands  which  are 
non-linear.  It  is  believed  that  the  theory  of  nonlinear  regression 
may  bo  brought  to  bear  on  this  problem. 


7.  COMPARISON  WITH  AN  ALTERNATIVE  METHOD 

In  previous  sections,  we  have  shown  how  to  obtain  a  cor.fidencr 
interval  for  E(W)  roughly  centered  about 
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where  N  is  the  number  of  busy  cycles  observed.  Me  know  by  the  strong 
law  of  large  numbers  that  as  N  -*■  «•  , 


N 


?  Y.(1> 


N  fa  k 


E{Y^^ }  -  E{a^}  E{W]  a.e. 


and 


f  I  V  E{ai} 

k^l 


a.e. 


so  that 


WN  +  E{W}  a.e. 
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Hence,  the  estimator  V?N  1b  consistent.  Given  any  e  >  0,  the 


Urn  PC | W  -  E(W) |  <  e  }  -  1. 

N-t«  N 


An  alternative  "stlraator  for  E{W)  is 


.  H-l 

%  -  a  l  w 

M  M  j-o  1 


For  many  queueing  syotems,  e.g.  the  M/G/l  and  GI/M/1  queues.  ia 

that  estimator  among  the  claaa  of  estimators 


.  M-l 

wj,  * w—  y  w. 

M,m  M-m  /—  j 

*  j-=m  J 


m  •=  0,1,2, .. .  ,M-1 


which  minimizes  the  "mean  square  error" 


E { (W '  -  E(W})4} 

M  ,m 


for  all  large  values  of  M,  see  BL0MQV1ST  (1970)  . 

Now  let  Nq  be  a  fixed  positive  integer,  and  let 


H0  -  Etoj) 


The  two  estimators  W  and  W',..  ,  require,  on  the  average,  the 

N0 

same  length  of  simulation,  and  it  is  of  interest  to  compare  their 

respective  differences  from  E{W}.  For  a  fixed  >  0,  define 
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P i  P{|WN  -  E{W) |  1  c0), 

P'  H  P{^[M0]  "  E{W)I  -  e0K 

By  replicating  a  number  of  simulation  experiments.  It  is  possible  to 

obtain  point  estimates  and  confidence  intervals  for  the  par.imeters 

P  and  P*.  To  illustrate,  400  Independent  simulation  replications 

were  made  for  the  M/M/1  example  of  Section  6,  with  X  5  and 

u  ■  10.  In  this  case  E{a^}  ■  2  and  we  set  Nq  ■  100  and  Mq  m  200. 

Point  estimates  and  confidence  Intervals  for  p  and  p'  are  shown 

in  Table  6  for  three  different  values  of  £q.  We  see  from  Table  6 

that  there  are  no  significant  differences  between  p  and  p'  , 

implying  that  one  can  generally  expect  comparable  "accuracy"  from 

the  estimators  WM  and  Wr'  , .  Similar  results  are  obtained  for 
N0  lM0J 

other  values  of  Nq  and  Mq  -  E{a^}  N„.  These  findings  lend 
respectability  to  the  methods  outlined  in  the  previous  sections,  for 
those  methods  would  be  of  questionable  value  if  they  were  based  on 
estimators  less  "accurate"  than  commonly  used  alternatives.  Of 
course  the  main  advantage  of  our  methods  over  alternatives  is  the 
ability  to  use  statistical  analysis  appropriate  tor  sequences  of 
i.i.d.  r.v.'s. 
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TABLE  6 

ESTIMATES  FOR  p  i  P{|«N  “EM  I  i.  E0  AND  p*  E  P{ I ”[mq]  "E  W  I  -  E0 
(400  rep’ications  with  X  -  5,  u  ■  10,  Nq  ■  100,  and  *  200;  E{W}  ■  0.1) 


Point  Estimates 

90%  Confidence  Intervals 

*•0 

P  P* 

P  P' 

0.01 

.220  .225 

L . 186,  .254]  [.191,  .259] 

0.02 

.455  .423 

[.414,  .496]  [.382,  .464] 

0.03 

.625  .623 

[.585,  .665]  [.583,  .663] 
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